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Abstract. Nonlinear oscillations of a 3D radial symmetric Bose-Einstein condensate 
under periodic variation in time of the atomic scattering length have been studied 
analytically and numerically. The time-dependent variational approach is used for 
the analysis of the characteristics of nonlinear resonances in the oscillations of the 
condensate. The bistability in oscillations of the BEC width is invistigated. The 
dependence of the BEC collapse threshold on the drive amplitude and parameters of 
the condensate and trap is found. Predictions of the theory are confirmed by numerical 
simulations of the full Gross-Pitaevski equation. 
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1. Introduction 

Investigation of collective excitations of a Bose-Einstein condensate (BEC) has attracted 
a great attention for last years pQ. Among the problems of these studies a dynamical 
response of the condensate to temporal perturbations in it as a whole is of natural 
interest. Many efforts have been devoted to the analysis of the effect of temporal 
modulations in a trap potential on oscillations of the condensate. The parametric and 
nonlinear resonances in BEC oscillations have been studied in works El HJ- The 
dynamics of the BEC under rapid strong perturbations of the trap potential has been 
considered in [SI IE], where the stabilization of the condensate is shown to exist. Since 
the dynamics is nonlinear the chaotic oscillations are also possible. Recently excitations 
of multimode dynamics and stochastization in two component condensate oscillations 
were observed experimentally in the work [Jj. The collapse induced by the noise in the 
trap has been studied in jS]. 

It should be noted that the effects of temporal modulations of a trap potential on 
the condensate and noncondensate parts of the cloud are similar. Interesting nontrivial 
dynamics of the condensate is also observed under temporal variation of the atomic 
scattering length by means of the Feschbach resonances [Hj- As distinct from the 
trap modulations, such variations affect on the condensed and thermal parts of an 
atomic cloud differently. The dynamics of the condensate is governed by the Gross- 
Pitaevskii equation with varying in time nonlinear term. The temporal modulations 
of a trap potential leads to the appearance of time dependent linear potential in the 
nonlinear Schrodinger equation. As to the modulation of the scattering length, a s , it 
leads to the Schrodinger equation with time- dependent nonlinear potential that opens 
a new interesting area in the matter waves theory. Here we can mention the Bosenova 
phenomenon occurring under rapid change of the scattering length value from positive 
(corresponding to the repulsion between atoms) to negative (corresponding to the 
attraction between atoms) [TUj . 

Periodic variation in time of the scattering length causes appearence of new 
effects. One of them is pattern formation in 2D BEC (TTj (see related ID problem 
in the context of the nonlinear optics [T2~]). Another interesting phenomenon is the 
existence of resonances in the macroscopic quantum tunnelling between two tunnel- 
coupled condensates (T3] . It should be yet mentioned a possible dynamical stabilization 
of 2D BEC under rapid periodic variation of the atomic scattering length and 
existence of Feshbach resonances. Collective oscillations of a ID Bose-Einstein gas 
under different external perturbations of the trap and the scattering lenght have been 
considered in the work [Tfij . 

Oscillations of 2D BEC have been studied analytically and numerically in the work 
[T7j . Numerical study of the 3D BEC oscillations has been given in ^H] • At the same 
time analytical consideration of oscillations of 3D BEC is still lacking. A free attractive 
3D BEC is known to collapse at any number of atoms. Existence of a trap leads to 
another situation. For a given trap potential a 3D BEC goes to collapse if the number 
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of atoms exceeds some critical value. This has been checked experimentally in work 

The purpose of this work is to investigate analytically and numerically nonlinear 
oscillations of the 3D radial symmetric BEC under periodic variations of the atomic 
scattering length. We will employ the time-dependent variational approach and derive 
a nonlinear ordinary differential equatin (ODE) for evolution of the condensate width. 
Predictions for nonlinear resonances and instability of a condensate we will compare 
with numerical simulations of the full time dependent Gross-Pitaevskii (GP) equation. 

The article is organized as follows. In section El we derive the variational model to 
describe a 3D Bose-Einstein condensate in a trap potential. In section El based on the 
Van-Der-Pole method necessary equations are derived to describe main and parametric 
resonances. In section 0] numerical simulations as well as their analysis are given and 
finally in section El the main results are summarized. 



2. Variational Model 

The total wave function for a 3D Bose-Einstein condensate in a trap potential V tr (r) 
may be described by the Gross-Pitaevskii equation: 

h 2 

ih% = + V tr (r)^ + g(t)\^\ 2 ^. (1) 

2m 

Here Vt r (r) = muj 2 r 2 /2 is a trap potential, the factor g(t) = 4irh 2 a s /m and a s is the 
atomic scattering length. We will suppose that the time-dependent scattering length a s 
is constant to leading order, with 

a s = a a(l + hsm(Q't)), (2) 

where a stands for the dimensionless scattering length. 

Such variation can be obtained using for example varying in time magnetic field. 
So the atomic scattering length can be expressed as 

asit) =aB { 1 + B A Bl (t) ) ' (3) 

where is the asymptotic value of a s and the magnetic field B\ = B 2 sinf2t. Then far 
from the resonance, we can approximate a s (t) = a so a(l + hsin(Qt)), with 

A \ , a B AB 2 



As to experimental values of the parameters B and A Feshbach resonances have been 
observed in Na at 853 and 907 G j20j, in 7 Li at 725 G 21J and in 85 Rb at 164 G with 
A = 11 G (221- Another way of making variations in the nonlinear term is to vary the 
tight confined potential parameters of the low dimensional system. For a 2D pancake 
BEC it may be the frequency of trap in the direction of tightening and the transverse 
trap frequency in the case of a ID cigar-shaped BEC. 
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Under the scaling r = tu, p = r /l ,l = sj^;, U = ^^Ana s lQ, the Gross-Pitaevskii 
equation can be written in the dimensionless form 
1 o 2 

iU T + -AU ~Y U + a ^ + hsm(QT))\U\ 2 U = 0, (4) 

where Q = Q'/u and a — 1, a = — 1 corresponds to attractive or repulsive interactions 
between atoms in the condensate respectively. 

For the radially symmetrical case the governing equation takes the form 

iU T + \u pp + -U p -^U + a(l + h sin({lT))\U\ 2 U = 0. (5) 

Corresponding Lagrangian density is determined by the following expression 

W V) = \{U*U T - UU%) - \\U P \ 2 - \p 2 \U\ 2 + |(1 + ^sin(C!r))|C/| 4 . (6) 

In deriving our variational model we proceed from the anzatz that the modulus of the 
total wave function, \U\, evolves in a self-similar way. So a trial function may be taken 
in the form 

ff ( ftT ) = A(r)fl(^)™ P (i^g + ^r)), (7) 

where the function R 2 f^yj describes initial BEC density distribution. 
To obtain equations for the wave packet parameters A(r), a(r), b(r), 0(r) we 
calculate the averaged Lagrangian L(r) = J d 3 pL (U(p, r), U*(p, r)) 
with the trial function (J7J). Its explicit expression takes the form 

(b T a-a T b + a 2 + b 2 )I 2 + 2(f) T n + ^-o-(l + hsm(nT))A 2 I 3 ) . (8) 



Hereafter the following designations for the constants n = J R(p) 2 p 2 dp, I 2 = 
J R(p) 2 p A dp, I 3 = / R(p) 4 p 2 dp and J 4 = / R p (p) 2 p 2 dp will be used. 

The Euler-Lagrange equations for the functional X(r) lead to the following 
equations for the scaling parameter a(r) and chirp b(r) : 

a T = b, 

h 3o-(l + hsm(ttT))A 2 a 3 I 3 

b T = - — a . 9 

J 2 a 3 2a 4 / 2 



Eliminating the parameter b from the above system and making substitution a = 
we get the following evolution equation for v : 



4 

v 



1 aPit) , s 

Vrr = -~V y, 10 

where 

P(t) = P + Pi sin(fir) , P = , ^1 = P h (11) 

and the quantity 

N 1 = — [ d 3 r\U\ 2 = A 2 v 3 [ r 2 drR 2 = A 2 v 3 n = (12) 
An J J In 
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determines the norm of the wave function U. 

It is easy to see that the above is a Hamiltonian system, with the scaling parameter 
v playing the role of a position variable, and the chirp b the conjugate momentum. Then 
this equation can be also rewritten as [3] 

d 



where 



dv V *ff( v )> ( 13 ) 
1 1 9 aP , , 

The dynamics of a condensate is different for the cases of positive and negative 
atomic scattering lengths. We will analyze these cases separately. First, for sake of 
completeness, let us consider P\ = 0, that is, the case when the scattering length is 
constant. 

2.1. Positive scattering length 

When cr = —1, i.e. a s0 > 0, the equilibrium points v of Eq.(fTU|) satisfy: 

vl = v + P . (15) 

This is a fifth-order algebraic polynomial problem that has only one positive real root. 
Its solution corresponds to a stable equilibrium point. When the interactions are strong 
Po 3> 1, one can neglect the term vq and obtain = Pq^ 5 . After substitution of 
v = t>o°^ + 8 into (US), we obtain 

V0 = 5(P )^-1- (16) 
This corresponds to the Thomas-Fermi regime, since we neglect the quantum pressure 
in comparison with parabolic potential and nonlinearity terms in the GP equation. 
The frequency of small oscillations near this point is 

^fi + dWl' ( 17 ) 



IjJr 

where it is assumed that P » 1 



V i0P n 4/5 



2.2. Negative scattering length 

A typical behavior of V e ff(v) for two cases, attractive and repulsive condensates, are 
presented in figure [TJ 

The equation for the equilibrium position 

^0 = ^3-^ PI) 

has either no, or two positive real solutions. The critical value of Po is given by 

P c = 5^4 = - 5350 - ( 18 ) 



Resonances in a trapped 3D Bose-Einstein condensate 



6 



When P Q > P c , there are no equilibrium points, and the condensate will collapse [23j. 
When Pq < P c there are two equilibrium points, one of them unstable (t>oi) and the 
other one stable (^02)- This means that collapse can be avoided when two conditions 
are satisfied: 

a) P < P c 

b) the initial condition: v 03 > v > t> 01 , where v 03 is defined by: V(v 03 ) = V(v Q i). 

In this region the motion should be periodic motion around the equilibrium point. 

When P <^ 1 for the equilibrium point of (jlUJI (the minimum of the potential V e ff, 
see figure^) we have the following approximation 

^0 = 1 - f ■ (19) 

In the interval < Pq < 0.5 the position of maximum of the potential function V e ff 
is well approximated by Vqi ~ P , the point v 03 ~ l/(v3P )- The frequency of small 
oscillations near v is 

(20) 

3. Periodic variation in scattering length 

One of the problems of interest in BEC is the study of the oscillation in the condensate 
subject to a periodic variation of the scattering length O H3 EE IH] • Such a variation, 
as was above mentioned, can be achieved experimentally by varying the magnetic field or 
using optically induced Feshbach resonances. In GP equation the temporal modulation 
of the scattering length of atoms corresponds to a time-dependent nonlinear coefficient. 

For the sinusoidal variations of the atomic scattering length the equation (|10j) 
takes the form: 

d 2 v 1 oPq aP l . 

-T^ + v = — 4 r sm{Qt). (21) 

We are looking for a solution of this equation in the form 




V = Vq + V-i, 

where V\ -C v q is a small deviation from the equilibrium position v . Substituting it into 
()21|) and expanding the terms of the equation near the equilibrium point v we get the 
equation for V\\ 

vi + uolvi + Av\ + Bv\ + ... = (C + D Vl + Ev\ + Fv\ + ...) sin(fit). (22) 

Above, we introduced the following constants: 
o~P 

ul = l + - 4—^, A = -Q/v 5 + 10aPo/vl B = 10/v 6 - 20aP /v 7 , 
v o 

C = -aPi/v%, D = AaPi/vl, E = -10<rPi/vg, F = 20aP l /v 7 . 
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In addition, we will assume that V\ and P\ are small quantities of the same order, say 
e. (Note that constants C, D, E and F depend on Pi.) Then, we can keep in equation 
(|22j) terms up to order e 3 : 

vt + ulv x + Av\ + Bv\ = {C + Dvt + Ev\) sin(ftt). (23) 

Now, we proceed on the examination of the resonance cases of equation (J2*3j) . We 
assume that 

UJq « -ft, 

q 

where p and q are mutually prime integers. Solution of equation (|2*H|) may be sought in 
the form of the expansion 

vt = a cos(ftt + 9) + hx {1) + h 2 x {2) + (24) 

where h is the amplitude of periodical perturbation (see equation (J2J)). Making use of 
the successive approximations method described in reference [21], in the case of main 
resonance (p = q — 1) we get the following set of equations for the amplitude, a and 
phase 9: 

^) cos W> ( 25 ) 



dt V 2ft 8ft 

dt = ~ n + 2^ smW + (sn " i2^n j a + m asm{6) - VB)> 

For the steady-state regime da/dt = d6/dt = 0. Then the amplitude frequency 
response can be easily found: 

(3a 3 ± 3Ea 2 ± AC + 8^70 = 0, (26) 

where 7 = cj — ft, | 7 |>C 1 and 

10A 2 

(3 = 3B- — r . dH) 
The signs (±) in (J26|) corresponds to different branches of the response (i.g. see figure 

EJ). 

System (|2~5|) can be transformed into hamiltonian form by changing the variable 
a = y/Q. In this case dQ/dt = dH/d9, and d9/dt = —dH/dQ, where 

H(Q, 9) = -^^Qsin(9) - ^Q 3 ' 2 sm(9) + (ft - u )Q - ^Q 2 . (27) 

In the case of parametric resonance (p = 1, q = 2, v 1 = a cos(|ftt + #)) the equations 
for a and 9 have the form: 

^ = -4Lcos(2fl) (28) 
dt 2ft v ; v ; 



«* " n« + ™ sm(20) + ( — - ^ ) a 2 (29) 



dt 2 2ft v y Uft 6^ft 
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Corresponding Hamiltonian is: 



f ) = sm(20) + (^-fi - co j Q - . (30) 

It follows that H is the integral of motion and there exists a separatrix H = 
separating finite trajectories corresponding to nonlinear resonances from infinite ones. 

Phase planes for two cases of attractive BEC (er = 1) with Pq = 0.2 and Pq = 0.4 
are shown in figure El As seen the portraits differ drastically. The reason is different 
signs of factor j3 entering equation ({25} for the cases P$ = 0.2 and Po = 0.4. 

Phase planes for two cases of repulsive BEC (er = —1) with P = 0.4 and P = 0.8 
are shown in figure El 

3.1. Dynamics of an attractive condensate under periodic a s 

Oscillations with the maximum amplitude for the case of the main resonance take place 
for separatrix H s « when 9 = n/2 (see equation (|27|l). The critical drive amplitude 
Pic providing the maximum amplitude of the attractive condensate width oscillations, 
a max may be obtained from (|26|). In order to obtain an expression for P lc let's rewrite 
above defined parameters C, D and E introducing new ones c, d, e as C = P±c, D = Pid 
and E = Pie. Then in the case of Q = u (7 = 0) we get 

I 3 I a 2 

Plc " |3ecw + ^r (31) 

"max 

The maximum amplitude a max is found as follows. As known the particle leaves the well 
(114)1 when its kinetic energy (vf t )/2 becomes comparable with the difference between 
the potential well bottom and separatrix, AV. This energy is: 

AV = V eff (v 01 )-V eff (v ), (32) 

where vq is stable and vqi is unstable equilibrium points (see subsection 12.2)) . For the 
case of slowly varying amplitude a and phase 9, we have from equation (j24j) : 

vu = —aQ sin(Qt + 9). (33) 

Equating critical value of the kinetical energy corresponding to the maximum amplitude 
of oscillations a max to AV we arrive to the following expression 



V2A~V 



(34) 



Another way of calculating the maximum amplitude, a max is to employ a direct 
expressiun for it 

a ma x = ^(v 03 - V i), (35) 

where the points Vq 3 and Vqi have been defined in 
subsection 12.21 
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In the region < Pq < 0.4 where the well is deep, we have the approximations 
Vox ~ Pq an d t>03 ~ -/p^ — j^Pq- Then in the indicated region of P f° r a-max we get a 
very good approximation 



V3P 



2V3 



Pn 



(36) 



4. The results of numerical simulations of the full Gross-Pitaevskii 
equation 

Based on the above-stated theory we can now proceed to numerical simulations to 
test calculations of the previous sections and carry out more detailed analysis. We 
discretize the problem in a standard way, with time step At and spatial step Ar, so Uj 
approximates U(jAr, kAr). More specifically we approximate the governing equation 
(JHJ) with the following second order accurate semi-implicit Crank-Nicholson scheme, 

i(U* +1 - U}) 1 



Ar 



2Ar 2 



2rjAr 



fc+1 
+1 



+ 



-j-(l + h8in(QT k ))a\Uj 



k\2 



(Uf + P fc+1 ) 



(37) 



The set of algebraic equations (|57j) is solved by the vectorial sweep method. 

In solving the full governing partial differential equation (PDE) (j5J) both for 
attractive BEC (er = +1) and repulsive one (a = — 1) we employ three different anzatzes 
as initial wave-packets : anzatz 1, called a Gaussian one fI5\ 12]; whith the trial 
function in (J7J) choosen as 



R 



exp 




and anzatzes 2 and 3 when the trial function R(-) — P(£) where R(£) is the solution of 
the equation 



ld 2 R ldR f n „ 3 .„ 

— z-7T H : — R + R - XR 

2 de e di 2 







(3f 



with A equal to —1 and +1 correspondingly. Such a choice enables us to use different 
initial wave packets. 

The parameters of the initial wave-packet are chosen to provide the given value of 
the norm N and P determined by Eqs.([T2|) and (fTTj). 

In PDE simulations the current values of the scaling parameter a(r) is calculated 
as follows : 



air) 



a , 



fd PP *\U(p,T) 



\ fdpp*\U(p,0\ 



(39) 
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Figure |U depicts profile of initial wave packets \U(p,r = 0)| corresponding to 
anzatzes 1, 2 and 3. Parameters of the wave packets are chosen to give the same 
value of the system parameter Pq = 0.2. 

Main wave packet parameters are given in Table below. 



wavepacket 


normN 


< r > 


anzatz 1 
anzatz 2 
anzatz 3 


0.01973 
0.02254 
0.01836 


1.06250 
1.07056 
1.06698 



As will be seen further in spite of some difference in the wave packets parameters 
their dynamics are atmost identical at small values of Pq. It should be also noted that 
anzatz 1 gives exact solution at Pq = 0, anzatz 2 does it at Pq = 0.372 and anzatz 3 
corresponds to a collapsing solution. 

4-1. Attractive condensate (a = +1) 

The results of numerical simulations of PDE and ODE models in the case of main 
resonance are presented below. 

Figure 03 depicts dynamics of the scaling parameter a of an attractive BEC under 
main resonance obtained from PDE and ODE calculations for two cases with: a) 
P = 0.2, Pi = 0.014 and b) P = 0.35, Pi = 0.00525. The frequency of the 
periodical perturbation is taken to be equal to the eigen frequency of the system ojq 
determined by 



In both cases of PDE simulations initial wave packet has gaussian form. The agreement 
between simulations of the full PDE and ODE in case a) is good enough for the time 
less than 100. In case b) when Pq = 0.35 there is clearly some discrepancy in the actual 
value of oscillation amplitude. It may be explained by the shift of the equilibrium point 
in full PDE calculations. 

In the case of main resonance the amplitude frequency response of the scaling 
parameter resonant oscillations is described by equation ()26)1 . The theoretical amplitude 
frequency response based on this equation and the one obtained from full PDE 
simulations for the case of Pq = 0.3 and Pi = 0.0045 are given in figure El In 
PDE calculations self similar anzatz (anzatz 2) from equation (ffi)]) is used as initial 
wave-packets. Each point of the plot results from PDE simulation of the scaling 
parameter evolution under main resonance starting from the stationary point Oq where 
the oscillations amplitude is equal to zero. 

The factor (3 in (J2T)|) plays an important role in the amplitude frequency response. 
Its sign determines the incline direction in the plot of the response and its value 
determines degree of the nonlinearity in the main resonance. Small values of the factor 
j3 means prevalence of the linear regime in the resonance. 




(40) 
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As was above shown collapse of the attractive BEC condensate can be avoided 
when P < P c (see ()18j0 . We have carried out full PDE simulations to study conditions 
when the main resonance makes the BEC to leave the potential well V e ff and the BEC 
begins to collapse. In terms of the effective Hamiltonian (}2Tj) this means that when in 
the resonance the amplitude of oscillations exceeds the value on the separatrix (H > 0), 
the condensate collapses. 

The dependence of the critical drive amplitude P\ c on P obtained from ODE 
(solid line) and PDE (filled squares) numerical simulations as well as the theoretical 
dependence ()31|) (dotted line) are shown in figure [7| The theoretical curve is obtained 
from equation (f3*T]) using approximation (jHHjl at < Po < 0.4 and expression (|H4"j) at 
0.4 > P < P\ c . As seen ODE simulations based on the equation (}2T]) is in a good 
agreement with PDE simulations of the full Gross-Pitaevski equation (JHJ) where self 
similar anzatz 2 (A = —1) is used as an initial wave packet. 

The dependences obtained from ODE and PDE simulations has a brightly expressed 
minimum at the point Po = 0.2. It may be explained by supposition that the resonance 
at this point is mainly linear. Indeed in accordance with equation (j26t ) the factor f3 
entering the folmula for the amplitude-frequency response goes to zero at Po ~ 0.277. It 
means that at this point the resonance is close to linear and small value of external drive 
Pi leads to large amplitudes of the BEC oscillations. Discrepancy between positions of 
the minimum in theoretical curve and the one obtained from numerical simulations may 
be explained by the fact that in deriving the governing equation ()22|) the terms of higher 
degrees than 3 are not taken in account. 

Next plot in figure |H1 depicts the dynamics of the parametric resonance (Q = 2u ) 
at the point Po = 0.36 when Pi = 0.051. Full PDE calculations were carried out with 
initial wave packet corresponding to anzatz 3 (A = 1). 

As seen at times t < 50 parametric resonance is proved to be suppressed in 
accordance with the prediction of the theory. 

4-2. Repulsive condensate (a = — 1) 

The dynamics of the scaling parameter a of a repulsive condensate under main resonance 
obtained from PDE and ODE calculations for the case Po = 0.4, Pi = 0.06 is shown in 
figure El The frequency of the periodical perturbation Q = ujq, where the eigen frequency 
u>o is determined as 



(for comparison see equation (J40j) ) . An initial wave packet has gaussian form (anzatz 
1). The agreement between simulations of the full PDE and ODE is excelent for the 
time less than 100. 

The dynamics of the repulsive BEC under parametric resonance resonance (Q = 
2ujq) obtained from full PDE calculations for Pq = 0.4 and Pi = 0.12 are shown in figure 




(41) 
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[TU1 Two cases are presented, the case when the initial wave packet : a) has self similar 
form (jHHj) with A = 1 (anzatz 3) and b) is gaussian (anzatz 3). 

5. Conclusion 

In this work we study the nonlinear oscillations of the 3D spherical symmetric Bose- 
Einstein condensate under periodic variation in time of the atomic scattering length. 
We employ the time-dependent variational approach and derive the equation for the 
dynamics of the condensate width. The equation has the form of the unharmonic 
oscillator under the singular periodic drive. Using this equation we analyze analytically 
the resonances characteristics for the cases of attractive and repulsive condensates 
. In particular we investigate the process of the attractive condensate collapse with 
N < N c under periodic variations of the scattering length. We find the possibility for 
the bistability in the oscillations condensate width depending on the detuning between 
drive frequency and the frequency of small oscillations of a condensate. 

We have shown that even when the number of the condensate atoms less than the 
critical, the condensate can collapse. Using the variational approach equation we derive 
the dependence of the driving force critical value P\ c versus initial parameter Pq. 

For the repulsive condensate case we study both main and parametric resonances 
in condensate oscillations. The full numerical simulations of the Gross-Pitaevskii 
equation with periodic nonlinear coefficient confirms predictions based on the variational 
approach. We find the deviations from predictions of the variational approach for the 
case when the number of atoms close to the critical one. 
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Figure 1. Effective potential used in describing the scaling parameter v. Two cases 
are presented : a) P = 0.4 (attractive BEC) and b) P = —0.4 (repulsive BEC). 



Figure 2. The phase planes for system {25J, <E3) at (a) P = 0.2 (Pi = 0.008, 
il = 1.936) and (b) P Q = 0.4 (Pi = 0.016, Q = 1.766) for a = 1 (attractive BEC), 
7 = and different H. 



Figure 3. The phase planes for the system |gS}, <G51) a * ( a ) Po = 0.4 (Pi = 0.06, 
fl = 2.066) and (b) P = 0.8 (Pi = 0.2, Q = 2.100) for a = -1, (repulsive BEC) 7 = 
and different H . 



Figure 4. Profile of initial wave packets \U(p, t = 0)| corresponding to anzatzes 1, 2 
and 3. Parameters of the wave packets are chosen to give the same value of the system 
parameter Pq — 0.2. 



Figure 5. Dynamics of the scaling parameter a of an attractive BEC under main 
resonance (f2 = u>o) obtained from PDE and ODE calculations, a) Po = 0.2, Pi = 0.014 
and b) P = 0.35, Pi = 0.00525. Solid line is for PDE and dotted line for ODE 
calculations. In both PDE simulations initial wavepackets are gaussian. 



Figure 6. Amplitude frequency responses obtained from full PDE simulations (dotted 
line) and from theory (solid line) for the case when Po = 0.372, Pi = 0.00372 and a = 1 
(attractive BEC). In PDE calculations anzatz 2 from equation 126|) is used as initial 
wave-packets. 



Figure 7. Dependence of the critical drive Pi c leading to collapse on Pq. Solid line 
is for ODE calculations and filled squares are for full PDE calculations. Self similar 
anzatz (anzatz 2) with A = —1 is chosen for initial wave packets in PDE simulations. 
Theoretical dependence do 1 ft is shown by dotted line. 



Figure 8. Dynamics of the scaling parameter a under parametric resonance (f2 = 2uo) 
obtained from full PDE calculations for Po = 0.36 and Pi = 0.051. As initial wave 
packet self similar anzatz 3 (A = 1) is used. 
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Figure 9. Dynamics of the scaling parameter a of a repulsive BEC under main 
resonance (fl = ujq) obtained from PDE and ODE calculations, Pq = 0.4, Pi = 0.06 
(ll> = 2.0664). Solid line is for PDE and dotted line for ODE calculations. An initial 
wavepacket is chosen to be gaussian. 



Figure 10. Dynamics of the repulsive BEC under parametric resonance (ft = 2w ) 
obtained from full PDE calculations when Po — 0.4 and Pi = 0.12. Two cases are 
presented : a ) when self similar anzatz 3 is used as an initial wave packet and b) when 
the initial wave packet is gaussian. 
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